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Abstract 

We propose and test a wavelet transform modulus maxima method for the au- 
tomated detection and extraction of coronal loops in extreme ultraviolet images of 
the solar corona. This method decomposes an image into a number of size scales and 
tracks enhanced power along each ridge corresponding to a coronal loop at each scale. 
We compare the results across scales and suggest the optimum set of parameters to 
maximise completeness while minimising detection of noise. For a test coronal image, 
we compare the global statistics (e.g., number of loops at each length) to previous 
automated coronal-loop detection algorithms. 

Keywords: Corona, Structures; Active Regions, Structure 



1. Introduction 

Historically, images of the Sun have always presented a myriad of features far in ad- 
vance of our physical understanding of their existence and evolution. New data with 
increased sensitivity, spatial- and temporal-resolution are continually challenging our 
theoretical models of the solar atmosphere. This general statement is especially 
true in the case of coronal loops as observed by the fleet of extreme ultraviolet 
(EUV) imagers launched since the mid 1990s. Beginning with the Extreme ultra- 
violet Imaging Telescope (EIT: Delaboudiniere et at, 1998), this list includes the 
Transition Region and Coronal Explorer (TRACE: Handy et al, 1999) and recently 
the Extreme Ultraviolet Imager (EUVI; Howard et al., 2008) onboard the twin Solar 
Terrestrial Earth Relations Observatory (STEREO: Kaiser et al, 2008) spacecraft. 
From these data, our current models suggest these coronal loops trace out hot coronal 
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plasma (wlMK) up to heights of around 50 Mm above the surface of the Sun, but 
questions remain over their temperature and density profile, temporal evolution, and 
3D structure. 

Studies of these open questions are limited by the ability to extract the loop 
system of interest from the hundreds of others typically present in a solar EUV 
image. Furthermore each individual loop as viewed by any current instrument is 
probably a collection of multiple strands, and there will be any number of these 
strands overlapping along the instrument line of sight. With the launch of STEREO, 
solar physicists can now view the corona from multiple angles, enabling the 3D 
reconstruction of coronal loops to try to overcome this problem. The largest problem 
in this reconstruction is the identification of the same feature from each spacecraft. In 
the ideal world, the same feature would be tracked in image sequences as observed 
from each spacecraft and the 3D reconstruction would become a mathematically 
simple problem. This would enable the scientist to proceed with studying the physical 
parameter of interest. Ideally this extraction should be automated in order to make 
the process instantly repeatable. From these issues, there arises a natural requirement 
for the automated detection of loops based on a statistical approach. Additionally, 
and perhaps most pressing, the expected data load from the Solar Dynamics Obser- 
vatory (SDO) will be overwhelming without automated feature recognition (McAteer 
et al, 2004, 2005b). 

Aschwanden et al. (2008) describe five such feature-recognition algorithms for 
coronal-loop identification and apply each algorithm to an example image from 
TRACE. This technique of comparing algorithms to each other, as well as to some 
ground truth, is the best method for rigorous testing of any feature recognition 
algorithm. In this paper we apply a sixth such technique to the same TRACE image 
as described in Section 2. The algorithm is described in detail in Section 3. The 
technique is called the 2D Wavelet-Transform Modulus Maxima (WTMM) method 
and was originally developed by Arneodo and colleagues as a multifractal analysis 
formalism (Arneodo, Decoster, and Roux, 2000; Arneodo et al, 2003; Khalil et at, 
2006) . It has been expanded to characterize the anisotropic nature of complex struc- 
tures (Khalil et al., 2006; Snow et al, 2008; Khalil et al, 2009) and also to perform 
an automated and objective segmentation of image features of interest from a noisy 
background (Khalil et al., 2007; Caddie et at, 2007). In solar physics, the technique 
has been used recently to study the complexity of solar active region magnetic fields 
(McAteer, Gallagher, and Ireland, 2005a; Conlon et al, 2008, 2010; Kestener et al., 
2010; McAteer, Gallagher, and Conlon, 2009), X-ray solar-flare emission (McAteer 
et al, 2007), and in tracking coronal mass ejections (Byrne et al., 2009). In Section 4 
we explain the natural advantages of applying this to coronal-loop identification and 
compare the global statistics against those identified in Aschwanden et al. (2008). 
Finally, in Section 5 we discuss the benefits and drawbacks of the WTMM algorithm 
and suggest future extensions. 

2. Observations 

In order to directly compare our results to those previously published, we have tested 
our algorithm on the same sub-image as in Aschwanden et al. (2008). This image was 
observed on 19 May 1998, 22:21:43 UT, with an exposure time of 23.172 seconds with 
a passband centered on 17lA. The data were de-biased, flat-fielded, and de-spiked 
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Figure 1. Edge detection of coronal loops for the underlying test image. The image is 800x600 
pixels, corresponding to 296 X 222 Mm. Maxima chains were obtained from the wavelet trans- 
form modulus maxima, using the analyzing wavelets defined in Equation (2) at a size scale of 
~ 11 pixels. 



to remove cosmic-ray spikes from the image. The resulting image is shown with our 
identified loops in Figure 1 and presents three of the main problems with detecting 
coronal loops: Firstly, they are curvi-linear structures: point sources, straight lines, 
and areas are well behaved in the [x, y] CCD plane; curves are much more difficult. 
Secondly, many of the loops rise almost vertically through the solar atmosphere. The 
rapid density drop off in the corona (with a scale height of about 50 Mm) results in 
a distinct lack of loop tops along many of these features. This makes it difficult to 
track from footpoint to footpoint. Thirdly there is a collection of small-scale features 
near the bunch of footpoints, which are typically not of interest in studies of coronal 
loop systems. 



3. Method 

Our method of coronal loop identification is based on the 2D Wavelet Transform 
Modulus Maxima method. The continuous nature of the wavelet transform allows 
us to scan all size scales continuously in order to take full advantage of the space- 
scale information available. This allows us to perform the segmentation of objects of 
interest in total objectivity, without any prior knowledge on the size or morphology 
of the objects. 

Image segmentation with continuous wavelets is based on the derivative of a 
2D smoothing function (filter) acting as an "edge detector". Let us consider two 
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wavelets that are, respectively, the partial derivatives with respect to x and y of a 
2D smoothing (Gaussian) function, 

^Gau(^)=e-^ + ^/ 2 = e-W 2 /2, (1) 

namely 

ihfay) = d(t>Ga,u(x,y)/dx 

and (2) 
i>2{x,y) = 94> Gaa {x,y)/dy. 

For any function f(x,y) 6 L 2 (EL) (where L 2 (R) consists of all square-integrable 
functions), the continuous wavelet transform of / with respect to ip\ and ip2 is 
expressed as a vector (Mallat and Zhong, 1992; Mallat and Hwang, 1992): 



TV[/](b,a) = 
' T in [/] = a~ 2 J d 2 x fa (a-^x - b))/(x) 
Tf 2 [/] = a~ 2 J d 2 x fa (a" 1 (x - b)) /(x) 
1 = V{T 0Ga J/](b,a)} = V{0 Gau , b , a * /}. 



(3) 



Thus, Equation (3) amounts to defining the 2D wavelet transform as the gradient 
vector of /(x) smoothed by dilated versions (^>Gau(a -1 x) of the Gaussian filter. The 
wavelet transform can be written in terms of its modulus, M,p[f}(h, a) and argument, 
A4/](b,a) 



M, 



> [/] (b, a) = ^/(T^[f}(b,a)) 2 + (T^[f](h,a)) 2 , (4) 



A^[f}(h,a) = Arg(T^[/](b, a ) + iT V)2 [/](b, a )). (5) 

The modulus maxima of the wavelet-transform, or intensity gradient maxima, are 
defined by the positions where the modulus of the wavelet transform, A4^,[f](h,a), 
i.e., the gradient, is locally maximal. These WTMM are automatically organized as 
maxima chains which act as contour lines of the smoothed image at the considered 
scales. This is a non-trivial process, which has been discussed and solved empirically 
by Arneodo, Decoster, and Roux (2000). At a given scale, the algorithm scans all the 
boundary lines that correspond to the highest values of the gradient, i.e., the maxima 
chains. For each size scale considered, the algorithm outputs the [x, y] pixel location 
chain of each edge detected. Each chain then corresponds to a single extracted edge 
in the image. We identify each edge as a single coronal loop, as shown in Figure 1. 
This process is repeated continuously at all scales. Futher detail on the actual code 
is given in the Appendix. 



4. Results 

In this section we describe the nature of the WTMM method to trace out and connect 
edges as one of the major advantages of this method. We compare automatically 
detected loops against those manually identified in Aschwanden et al. (2008). In 
order to quantitatively compare this technique against those previously published we 
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Figure 2. A complete coronal loop, (a) Thefrr, y] location of each pixel along the loop from 
one footpoint. The colour scheme traces the loop from the first footpoint (green), through the 
loop top (red) to the second footpoint (yellow). The axis are in pixel units, where one pixel 
~ 370 km. (b) The Modulus [Equation (4)] and Angle [Equation (5)] information at each 
point along the loop with the same colour scheme as in (a), where a positive (negative) angle 
corresponds to a counter-clockwise (clockwise) rotation from (along the positive ^-direction) 




Figure 3. An incomplete coronal loop, (a) The [x, y] location of each pixel along the loop 
from one footpoint. The colour scheme traces the loop from the first footpoint (yellow), up 
through the atmosphere (red, green) where it presumably connects to another footpoint, pos- 
sibly elsewhere in the image. The axis are in pixel units, where one pixel 370 km. (b) The 
Modulus [Equation (4)] and Angle [Equation (5)] information at each point along the loop 
with the same colour scheme as in (a), where a positive (negative) angle corresponds to a 
counter-clockwise (clockwise) rotation from (along the positive ^-direction) 



also collate a number of global image statistics: the total number of loops detected 
is an indication of the completeness of the algorithm: the maximum loop length is 
an indication of the upper limit on the detection: the scaling index of the cumulative 
distribution shows how the algorithm performs across size scales. 



4.1. Loop Tracing 



Figure 2a shows an individual loop corresponding to one of the large complete coronal 
loops near the center of the image. The modulus and angle information used to detect 
this loop is shown in Figure 2b with the same color scale as Figure 2a. As the loop 
is traced from one footpoint to the other, both the modulus and angle information 
change smoothly. The modulus is strongest at the footpoints (where the intensity in 
the original image is strongest) and drops off by a factor of five at the loop tops. The 
angle starts off at 0° at the right-most footpoint (i.e., pointing to the right of the 
image, along the £-axis), increases to tt/2 at the loop top (i.e., pointing straight up 
along the y -axis), and increases to it at the second footprint (i.e., pointing to the 



SQLA1029R3_McAteer.tex; 17/02/2010; 1:34; p. 5 



Solar Physics 



Loops Scale S Scale 6 Scale 12 



□□B 

Figure 4. Comparison of (left to right) the smooth-subtracted data, manually identified loops, 
edges from scale 2 of the WTMM, edges from scale 6 of WTMM, edges from scale 12 of 
WTMM. Top row is for a subimage extracted at the leftmost footpoints of the loop system 
consisting of mostly vertical loops, bottom row for a subimage extracted near the rightmost 
footpoints of the system and consists of more curved loops. Each sub-image is 100 X 150 pixels 
(37 000 x 55 500 km) 

left, along the x -axis). It is not a completely smooth transition; this is most evident 
at the loop top, where the angle varies more rapidly. This variation of angle and 
modulus hence provides two parameters for loop identification. Most importantly, 
although the modulus may drop off dramatically at the loop top (a manifestation of 
the decreased signal), the variation in angle is much smoother. A second example, 
for tracing a partial loop, is displayed in Figure 3 for a loop leg in the lower right 
of the image. In this case, the angle stays constant (w n/2) along the leg, while the 
modulus drops off with increasing distance from the footpoint. 

4.2. Comparison of Loop Locations 

A more detailed analysis of extracted loops can be carried out by comparing the 
edges against those loop manually identified in Aschwanden et al. (2008). Figure 4 
shows such a comparison for two subimages of the data. These regions are picked as 
they illustrate two important spatial regimes of loop identification and demonstrate 
some of the main positive aspects and drawbacks of the WTMM method. Figure 4 
(top) contains a bunch of loop legs. These are mostly linear and display a relatively 
constant intensity along each loop inside this window. These two aspects should 
make these features easier to detect however none of the WTMM scales completely 
extract the manual coordinates. Each of the scales displayed pick out seven -eight of 
thirteen manually detected loops. The largest scale displayed (scale 12; right column) 
is unable to pick out any of the weaker, shorter, loops however does produce smooth 
edges. The smallest scale (scale 2; middle column) identifies more of the smaller 
loops, but at the expense of producing a more ragged edge. Scale 6 still contains 
many of small loops, but also produces a smooth edge. Scale 6 is also the only scale 
which captures the loops in the bottom right of this subimage. Figure 4 (bottom) 
shows a selection of the curvi-linear features. As these are curves in [x, y] pixel space 
and show a drop in intensity, they are normally more difficult to extract. However 
the ability of the WTMM method to chain in angle space assists in overcoming this 
problem. Of the seven curved, manually-detected features, each scale reproduces 
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Figure 5. The completeness and upper limits across each wavelet scale. The bottom axis 
corresponds to the wavelet scale order; the top axis is the same value converted to pixel units. 
The total number of loops with size greater than 70 pixels [TV], is plotted as a solid line and 
asterisks according to the left axis. The solid horizontal line is plotted at the desired complete 
value of N = 154. The maximum loop length [L] is plotted as a dashed line with asterisks 
according to the right axis. The dashed horizontal line is plotted at the desired upper limit 
of L = 463. Both desired values are from the manually identified loops in Aschwanden et al. 
(2008) 



four -five. The large, but weak, loop near the center of the subimage is not detected 
at any scale. The wavelet transform power lies below the threshold for idenfication. 
Lowering our threshold to extract this loop results in pulling out an unmanageable 
number of weak non-loop features. Scale 2 and scale 6 are similar; the main difference 
is the ability of scale 6 to begin to detect the loop in the top left. At scale 12, the 
smoothing window is so large that the tracking algorithm jumps from one large loop 
to a second large loop. Clearly there is a trade off to be made between completeness 
and quality of extracted edges. 

4.3. Number of Loops and Maximum Loop Length 

A full quantitative comparison of the positive and negatives of each scale is best 
achieved by comparing the global statistics. The total number of detected loops and 
maximum loop length at each wavelet scale size is displayed in Figure 5. The bottom 
axis shows the wavelet scale, with the corresponding scale in pixels displayed on the 
top axis. The solid line shows how the total number of loops greater than 70 pixels 
(left axis) varies with each wavelet scale size. The horizontal solid line shows the 
ground truth (i.e., manually extracted) of N = 154. It is clear that the wavelet 
algorithm consistently underperforms, but it is noticeable that the plot peaks at 
scale size 6. At this scale we can reach a completeness ratio of 0.91. The dashed 
line shows how the maximum loop length (right axis) varies with each scale size. 
The horizontal dashed line shows the ground truth (i.e., manually extracted) value 
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Figure 6. The cumulative distribution of loops with length greater than 70 pixels 

of L = 463. As expected, the ability of the algorithm to detect the largest loops 
increases with increased wavelet scale size. It underperforms at small scale size, and 
over-detects larger structures at large scales sizes. The best performance is also at 
wavelet scale 6, where we reach a detection ratio of 1.05. It is comforting to note 
that we achieve maximum accuracy and completeness in these two parameters at 
the same scale size. 

4.4. Scaling Index of Cumulative Frequency Distribution 

The cumulative distribution of loops with length greater than 70 pixels is displayed 
in Figure 6 for scale 6 (« 11 pixels). The scaling index (/3) of this plot is indicative 
of the distribution of loop lengths in the image and reflects the tendency of smaller 
loops of naturally outnumber larger loops. Aschwanden et al. (2008) report a ground 
truth of P = —2.8 for the manually traced loops, with the five algorithms producing 
values between -2.0 and -3.2. For the WTMM algorithm in Figure 6, /3 = -2.78±0.1. 

5. Conclusions and Future Work 

EUV images of the corona provide a multitude of information regarding the plasma 
properties of coronal loop systems. A number of algorithms currently exist that 
attempt to extract the locations of these loops automatically. These algorithms are 
discussed in detail in Aschwanden et al. (2008) but here we compare a brief outline 
of each one against our WTMM method. 

The Oriented Connectivity Method (OCM: Lee, Newman, and Gary, 2006a) con- 
sists of a preprocessing step to remove non-loop candidates, a linkage step based 
on magnetic-field extrapolation guide, and post-processing to spline fit and link 
loop segments. The Dynamic Aperture-Based Loop Segmentation Method (DAM: 
Lee, Newman, and Gary, 2006b ) replaces the linkage step of the OCM with a 
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search for loops by connecting pixels which have a similar Gaussian cross sectional 
profile (Carcdeo et at, 2003) and orientation. The Oriented-Directivity Loop Tracing 
Method (ODM: Aschwanden et al, 2008 ) is essentially a local-directivity version of 
the OCM, tracing loops locally in two directions from their tops to their footpoints. 
The Ridge Detection by Automated Scaling (RAS: fnhester, Feng, and Wiegelmann, 
2008 ) is is a multiscale ridgel extension of the OCM. Finally, the Unbiased Detection 
of Curvi-Linear Structures (UDM: Steger, 1996 ) is a more generic method consisting 
the determination of the centroid of the loop structures from a second derivative 
perpendicular to the loop, and extended (Raghupathy, 2004) to connect structrues 
using a generalised Radon transform. Aschwanden et al. (2008) show that the OCM 
and DAM both successfully extract the large loop features. The ODM, RAS, and 
UDM codes contain many more free parameters and extracted more segmented small 
loop-like features. 

The WTMM method presented here offers many natural advantages over other 
techniques: Firstly, it has naturally directional linkage. The algorithm works by 
tracking along edges (perpendicular to the gradient) hence the resulting edges require 
minimum post-processing linkage. Secondly it is naturally multiscale, enabling a 
scientist to extract features in the size-scale range of interest. Essentially we negate 
the need for preprocessing by simply studying the scale at which the feature of 
interest occurs. These two features combine to provide a good localisation and a 
single response to each edge in the image. These advantages are evident in the good 
comparison of global statistics against existing algorithms from Aschwanden et al. 
(2008). We reach a completeness of 0.91 and detection ratio of 1.05, with a scaling 
index for the cumulative distribution of /3 = —2.78 ± 0.1. These three statistics 
compare favourably to those in Aschwanden et al. (2008). The WTMM contains a 
few free parameters which we have attempted to optimise for loop detection. The 
most important free parameter is the choice of the form of the mother wavelet. We 
choose the Derivative of Gaussian (DOG) as it is well studied mathematically and 
has previously been successful in tracking edges. We note that other, more naturally 
curvi-linear multiscale algorithms exist (e.g., curvelets, ridgelets) which may assist 
in providing better identification of individual coronal loops. The other main free 
parameters are the choice of thresholds in tracking the edges (see Appendix for 
details). Finally we propose that our sixth wavelet scale (corresponding to w 11 
pixels) seems to optimise our ability to pick out the small scale features, retain 
smoothness in the extracted coordinates, and best agrees with the expected global 
statistics 

An obvious extension of this work is to use the multiscale aspect in a soft- 
thresholding sense: currently we decide on an optiumum scale, there is the possiblity 
to instead decide on a best scaling index. This would consist of attempting to track 
each feature across scale. A noisy feature is expected to contain a lot of power at 
small scales, with very little power at larger scales. For a real feature, particularly 
the large complete coronal loops, there is expected to be much less change across 
scales. This may assist in removing non-loop features near the footpoints of the loops 
system. We also note our extracted loops, even at scale 6, are probably too ragged 
for stereoscopy and a degree of post-processing smoothing may still be neccessary. 
When applied to EUVI images from both STEREO ahead and behind data, this 
may allow for a 3D reconstruction of these loop systems with less manual labour. 
We expect algorithms such as the one studied here will be vital in the SDO era for 
the tracking of loops across a sequence of images, and hence in the expanding field 
of coronal seismology. 
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Appendix 

A. Light Weight Notations for WTMM Formulae 

Let f(x, y) be the input image to be analized and T the wavelet transform vector of 
/ at scale a. The wavelet transform components are 

fx = d x (f*<f> a ) (6) 
fy = d y (f*<j, a ) (7) 

The square modulus of the wavelet transform vector is M 2 = f 2 + f 2 . Let us note 
higher order derivative of / this way using multiple x- or ^-indices: d x f x = f xx , 

Oyfx fxy? &tc... 



B. WTMM Edge Definition 

The WTMM are defined as the locations of the points where there is a maximum of 
the wavelet transform modulus along the direction of the wavelet transform vector, 
i.e., WTMM are the location of the greatest slope in the f * (f>a landscape. The 
steepest slope line is not exactly orthogonal to the WTMM edge (only in particular 
cases). To find those points, we need to evaluate the scalar quantity (TV) defined as 
the dot product 

TV = V(Tl/ 2 ) • T (8) 

At each WTMM location, TV is zero and TV changes its sign when moving along 
the direction T and crossing the WTMM. This is not enough to clearly identify a 
maximum, we also require that the second derivative along the direction T must be 
strictly negative, i.e., 

N' = d x Nf x + d y Nf y < (9) 
The quantities TV and TV' are computed exactly: 

TV = V(M 2 ) . T = 2f 2 f xx + if x f y f xy + 2f 2 J yy . (10) 

TV' is given by: 

N' = 8 X Nf X + dyNfy = Afj 2 xx + Af y f 2 y +2f x f XXX + 2f y fyyy 

+4f 2 f 2 +4f 2 f 2 

1 J x J xy 1 J y J xy 
~^~&fxfyfxyfyy H~ &fxfyfxyfxx 
+Gf x fyf„y+6f,fy fxyy 



C. WTMM Computation Algorithm 

Finally, the full algorithm is summed up. It is essentially a FOR loop over pixel 
location. Note that the final edge image is made of the pixels containing a WTMM. 
The modulus at the "exact" WTMM location is adjusted using a polynomial fit using 
points along direction T. 
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Algorithm 1 WTMM edges computation algorithm 
Require: f(i,j) the input image 
Require: a scale parameter 

compute WT of / at scale a : T = V(/ * G„) 

compute iV and N' 

for pixel € image range do 

if N' < and N changes sign in 3 x 3-neighbourhood then 
pixel is labelled as a WTMM 

perform a third order polynomial interpolation along direction T to get an 
accurate wavelet transform modulus value at maximum 
end if 
end for 

WTMM edge image 



Several free parameters have been statistically and/or empirically adjusted over 
the years. The following have been used for this study: 

- Spatial resolution: Finite-size and edge effects. The minimum scale at which the 
wavelet transform is still resolved is seven pixels. On the other hand, in order to avoid 
artificial effects from the edges of the images, only the central 72% of the original 
wavelet transformed image should be kept for analysis. A methodical calculation of 
this parameter is carried out in Arneodo, Decoster, and Roux (2000) 

- Scale resolution. For an image of 1024 x 1024 pixels, usually 50 wavelet scales are 
considered, from seven pixels to ~ 200 pixels, in log 2 steps (i.e., very high resolution 
for small scales and low resolution at high scales). The seven-pixel lower limit restricts 
our resolution and is adopted to reduce the computational overhead. 
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